*Main data pathway globals, set directory to clean data folder
global results "S:\Project\DemoSos2\common\felles\JR_RG\DrVA\ResultsRev2\"
global cleandata "S:\Project\DemoSos2\common\felles\JR_RG\DrVA\CleanData\"
cd "$cleandata"

cap log close
clear all
set more off
set matsize 11000
set maxvar 120000


use  clean_patientlevel_file_any_rev3.dta, clear
keep if yr_swap >= 2005 & yr_swap<=2014
keep if swap_GPIDnew!=. 

* largest mobility group 
a2group, individual(swap_GPIDnew) unit(swap_GPIDprev) groupvar(pair)
bys pair: ge size_pair = _N
tab pair
drop if size_pair < 10 

keep if swap_GPIDnew >=55

ge death2 = mortality_2year
for num 25 35 45 55: ge death2_X =mortality_2year if swap_age>=X

label var death2_55 	"2-years mortality for 55+"
global index "death2_55"

/* estimate the FE */
global controls "i.yr_swap i.swap_age male"


foreach var in   $index {
	reghdfe `var' $controls, a(_fe_`var'=swap_GPIDnew fe_old_`var'=swap_GPIDprev) res(res_`var')  groupvar(mobil_`var')

	predict x_`var' if e(sample), xb
	ge fe2_`var' = _fe_`var' + res_`var' if e(sample) 
	bys swap_GPIDnew: egen fe_`var' = mean(fe2_`var') 
	ge n`var' = e(N) if e(sample)	
	drop _fe_`var'
}



*---------------------------------------------------------*
* OUTPUT INFO
*---------------------------------------------------------*
bys swap_GPIDnew: gen count_doc = _n

foreach x in $index  {
	bys swap_GPIDnew: egen var_`x' = var(fe2_`x') 					// variance of fixed effect within doc
	egen var_across_`x' = var(fe_`x') if count_doc==1 				// var of all average fixed effects (across doc)
	bys swap_GPIDnew:  egen count2_`x' = count(fe_`x') 				// count number of fams per new doc
	bys swap_GPIDnew:  egen count3_`x' = mean(count2_`x')
	egen shrink_a_`x' = mean(var_`x'/count3_`x') if count_doc==1 	//  this is the variance of each doc's FE divided by the number of indiv a doc sees
	g shrink_b_`x' = var_across_`x' - (var_`x'/count2_`x')
	g shrinkage_`x'  = shrink_b_`x' / (shrink_b_`x' + shrink_a_`x') 
	egen shrink_coef_`x' = mean(shrinkage_`x') 						//  this is the shrinkage coefficient
	g afe_`x' = (sqrt(shrink_coef_`x'))*fe_`x'

}


/* TABLES */
keep if count_doc==1
global fe "fe_death2_55"
label var fe_death2_55 		"2-years mortality for 55+"

eststo s1: estpost tabstat $fe, s(n p10 p25 p50 p75 p90 mean sd) columns(statistics)  
esttab s1 using "$results/fe_unadjusted_any.tex", replace fragment ///
	cells("count(fmt(0)) p10(fmt(3)) p25(fmt(3)) p50(fmt(3)) p75(fmt(3)) p90(fmt(3)) mean(fmt(3)) sd(fmt(3))") noobs nomtitle ///
	varlabels (`e(labels)')  varwidth(20) booktabs label

keep swap_GPIDnew fe_*
rename fe_death2_55 anyswap_fe
duplicates drop
save "$results/fe_addon_any.dta", replace
append using "$results/fe_addon2.dta"

twoway 	(kdensity anyswap_fe, lcolor(red) lpattern(dash) lwidth(thick) graphregion(color(white)) bgcolor(white) xtitle(Value Added)) ///
		(kdensity fe_death2_55, lcolor(blue) lpattern(solid) lwidth(thick)), legend(order(1 "VA - any swap" 2 "VA - exogenous swap")) 
		graph export  "$results\figure_a8.pdf",  replace

		

*---------------------------------------------------------*
******FE analisys****
*---------------------------------------------------------*
use  clean_patientlevel_file_any_rev3.dta, clear
keep if yr_swap >= 2005 & yr_swap<=2014
keep if swap_GPIDnew!=. 

drop antbarn bigcity birthorder
merge 1:1 lopenr using "FinalSet4", keepusing(birthorder) keep(matched)  nogen
merge 1:1 lopenr using kids_rev, keep(matched)  nogen

* age of death
drop yob age_death
ge yob = yr_swap -swap_age
ge age_death = yr_death-yob

 
* largest mobility group 
a2group, individual(swap_GPIDnew) unit(swap_GPIDprev) groupvar(pair)
bys pair: ge size_pair = _N
tab pair
drop if size_pair < 10 

keep if swap_GPIDnew >=55


merge m:1 swap_GPIDnew using "$results/fe_addon_any.dta"
drop _merge

ge death5 = mortality_5year
ge death2 = mortality_2year
for num 25 35 45 55: ge death5_X = mortality_5year if swap_age>=X
for num 25 35 45 55: ge death2_X = mortality_2year if swap_age>=X
for var heart cancer ext resp mental accident suicide homicide: ge X2y_55 = X_2y if swap_age>=55


* 1 - observables

label var death2_55 	"2-years mortality for 55+"

label var heart2y_55 	"Mortality: heart conditions"
label var cancer2y_55 	"Mortality: cancer"
label var ext2y_55 		"Mortality: external conditions"
label var resp2y_55 	"Mortality: respiratory conditions"
label var mental2y_55 	"Mortality: mental health conditions"
label var accident2y_55	"Mortality: accident conditions"

bys swap_GPIDnew: ge n = _n
* recode main VA measure (more VA - better doc, less mortality) 
ge newva = -anyswap_fe
local va "newva" 
su `va', d
gen stdva = (`va'-r(mean))/r(sd)
local va "stdva"
lab var stdva "Standardized VA"

keep if swap_GPIDnew!=. & stdva!=. & swap_age>= 55

egen any_chronic_any0 = rmax(dcancer_any ddigest_chronic_any dheart_chronic_any dendoc_chronic_any dany_cardiogp_any) if yr_swap>=2008
egen chronic_any0 = rsum(cancer_any digest_chronic_any heart_chronic_any endoc_chronic_any any_cardiogp_any) if yr_swap>=2008 
drop hs0 college0
gen hs0 = inlist(bu_1_any,1,2) if  bu_1_any!=. & bu_1_any!=9
gen college0 = inlist(bu_1_any,6,7,8) if bu_1_any!=. & bu_1_any!=9

local baseline "college0 married_any NORborn antbarn tot_income_any dwoverfor_any total_days_sl_any total_any any_chronic_any0 chronic_any0 reimburse_any dnpr_any npr_any "

global controls "i.yr_swap i.swap_age male" 

*Run balance test with previous GP FE and controls, standardized, append
foreach var in `baseline' {
	eststo `var': qui areg `var' stdva $controls,  a(swap_GPIDprev) cl(swap_GPIDnew) 
	su `var' if e(sample)
	estadd scalar mean = r(mean)
	eststo `var'
}

esttab `baseline' using "$results/table_a8.tex", ///
	booktabs legend replace keep(stdva)     ///
	label star(* 0.1 ** 0.05 *** 0.01) stats(emp mean N, fmt(%9.3g)) b(3) se(3) ///
	mtitles("College Degree" "Married" "Born in Norway" "Nb. Childrens" "Income" "Any Welfare Benefits" "Sick Leave" "Nb. GP Visits" "Any Chron Cond" "Chron Cond" "Reimbursements" "Any Hospital." "Nb. Hospital.")	
eststo clear


* doc characteristics
ge dr_long_list = dr_list_length_str_exog_swap > 1100 if dr_list_length_str_exog_swap!=.
ge dr_old = dr_age_str_exog_swap >=40 if dr_age_str_exog_swap!=.
ge dr_max_list =  dr_max_str_exog_swap> 1200 if  dr_max_str_exog_swap!=.

replace  sl_days5  = 0 if  sl_days5 == . & tot_income5!=.

global dr_any "dr_age_swap dr_spec_swap dr_list_length_swap dr_male_swap dr_share_swap dr_group_swap dr_max_swap swap_nb_ref_short anyswap_fe"

global patient_any "total_any dlabtest_any dprocedure_any dsl_any"	



*-------------------------------------------------------------------
* characteristics of gps
*-------------------------------------------------------------------
preserve
keep if swap_age >= 55 & stdva!=.
collapse $dr_any $patient_any, by(swap_GPIDnew)
for var $patient_any: replace X = 0 if X==.
eststo s3: estpost tabstat  $dr_any $patient_any if anyswap_fe!=., s(mean sd) columns(statistics)
count if anyswap_fe!=.
restore 

lab var	total_any "GP consultations"
lab var dlabtest_any 	"Lab tests"
lab var dprocedure_any 	"Surgery"
lab var dsl_any  		"Sick Leave"

lab var dr_age_swap 	"GP Age" 
lab var dr_spec_swap 	"GP has specialization" 
lab var dr_list_length_swap "GP list lenght"
lab var dr_male_swap  	"GP is male"
lab var dr_share_swap 	"GP has a shared list"
lab var dr_group_swap 	"GP works in shared office"
lab var dr_max_swap   	"GP max list lenght"
lab var anyswap_fe		"VA: 2 years mortality rate 55+"
  
esttab s3 using "$results/tablea3_col5_6.tex", replace fragment cells("mean(fmt(2)) sd(fmt(2))") noobs nomtitle ///
varlabels (`e(labels)') varwidth(20) booktabs label mgroups("Exogenous Swaps" "Any Swap", pattern(1 1)) collabels("Mean" "SD")

*-------------------------------------------------------------------
* characteristics of patients
*-------------------------------------------------------------------
drop *2001* *2002* *2003* *2004*

ge hs = edu_yrs <= 13
ge college = edu_yrs >= 16

for var heart2y_55 cancer2y_55 ext2y_55 resp2y_55: ge _X = X if mortality_2year == 1

global health 	"death2_55  _heart2y_55 _cancer2y_55 _ext2y_55 _resp2y_55"
global demo 	"swap_age male NORborn hs edu_yrs married_any"
global use		"total_any dlabtest_any dprocedure_any dsl_any dreimburse_any dnpr_any"
global inc	 	"lab_income_any tot_income_any dwoverfor_any total_days_sl_any"

replace total_days_sl_any = 0 if  total_days_sl_any==.

* characteristics of indiv (all) 
estpost tabstat $health $demo $use $inc if  str_exog_age >= 55, s(mean sd) columns(statistics)  
mat A =[e(mean)',e(sd)']

count 
esttab matrix(A, fmt(3)) using "$results/tablea5_col5_6.tex", mlabel("\multicolumn{2}{c}{Any}\\") replace 




	